function triangleGrid=createTriangleGrid(sites,triangles)

triangleGrid.minX=min(sites(:,1));
triangleGrid.maxX=max(sites(:,1));
triangleGrid.minY=min(sites(:,2));
triangleGrid.maxY=max(sites(:,2));
numTriangles=size(triangles,1);

triangleGrid.intervalX=(triangleGrid.maxX-triangleGrid.minX)/sqrt(numTriangles);
triangleGrid.intervalY=(triangleGrid.maxY-triangleGrid.minY)/sqrt(numTriangles);

triangleGrid.numX=round((triangleGrid.maxX-triangleGrid.minX)/triangleGrid.intervalX);
triangleGrid.numY=round((triangleGrid.maxY-triangleGrid.minY)/triangleGrid.intervalY);

% siteGrid=cell(triangleGrid.numX, triangleGrid.numY);
% numSites=size(sites,1);
% for s=1:numSites
%     [x_idx, y_idx]=findIdx(sites(s,1), sites(s,2), triangleGrid);
%     siteGrid{x_idx, y_idx}=[siteGrid{x_idx, y_idx}, s];
% end

triangleGrid.grid=cell(triangleGrid.numX, triangleGrid.numY);
numTriangles=size(triangles,1);
for t=1:numTriangles
    [x_idx, y_idx]=findIdx(sites(triangles(t,:)',1), sites(triangles(t,:)',2), triangleGrid);
    [BBox_x, BBox_y]=meshgrid(min(x_idx):1:max(x_idx),min(y_idx):1:max(y_idx));
    theLength=size(BBox_x,1)*size(BBox_x,2);
    for l=1:theLength
        triangleGrid.grid{BBox_x(l), BBox_y(l)}=[triangleGrid.grid{BBox_x(l), BBox_y(l)},t];
    end
end



